/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.lucene.spatial3d;

import com.carrotsearch.randomizedtesting.generators.RandomNumbers;
import com.carrotsearch.randomizedtesting.generators.RandomPicks;
import java.io.IOException;
import java.io.PrintWriter;
import java.io.StringWriter;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashSet;
import java.util.List;
import java.util.Random;
import java.util.Set;
import java.util.function.IntSupplier;
import org.apache.lucene.codecs.Codec;
import org.apache.lucene.codecs.FilterCodec;
import org.apache.lucene.codecs.PointsFormat;
import org.apache.lucene.codecs.PointsReader;
import org.apache.lucene.codecs.PointsWriter;
import org.apache.lucene.codecs.lucene90.Lucene90PointsReader;
import org.apache.lucene.codecs.lucene90.Lucene90PointsWriter;
import org.apache.lucene.document.Document;
import org.apache.lucene.document.Field;
import org.apache.lucene.document.NumericDocValuesField;
import org.apache.lucene.document.StringField;
import org.apache.lucene.geo.Polygon;
import org.apache.lucene.geo.Rectangle;
import org.apache.lucene.index.DirectoryReader;
import org.apache.lucene.index.IndexReader;
import org.apache.lucene.index.IndexWriter;
import org.apache.lucene.index.IndexWriterConfig;
import org.apache.lucene.index.LeafReader;
import org.apache.lucene.index.MultiDocValues;
import org.apache.lucene.index.NumericDocValues;
import org.apache.lucene.index.PointValues.IntersectVisitor;
import org.apache.lucene.index.PointValues.Relation;
import org.apache.lucene.index.ReaderUtil;
import org.apache.lucene.index.SegmentReadState;
import org.apache.lucene.index.SegmentWriteState;
import org.apache.lucene.index.Term;
import org.apache.lucene.search.IndexSearcher;
import org.apache.lucene.search.Query;
import org.apache.lucene.spatial3d.geom.GeoArea;
import org.apache.lucene.spatial3d.geom.GeoAreaFactory;
import org.apache.lucene.spatial3d.geom.GeoBBoxFactory;
import org.apache.lucene.spatial3d.geom.GeoCircleFactory;
import org.apache.lucene.spatial3d.geom.GeoPathFactory;
import org.apache.lucene.spatial3d.geom.GeoPoint;
import org.apache.lucene.spatial3d.geom.GeoPolygon;
import org.apache.lucene.spatial3d.geom.GeoPolygonFactory;
import org.apache.lucene.spatial3d.geom.GeoShape;
import org.apache.lucene.spatial3d.geom.Plane;
import org.apache.lucene.spatial3d.geom.PlanetModel;
import org.apache.lucene.spatial3d.geom.SidedPlane;
import org.apache.lucene.spatial3d.geom.XYZBounds;
import org.apache.lucene.spatial3d.geom.XYZSolid;
import org.apache.lucene.spatial3d.geom.XYZSolidFactory;
import org.apache.lucene.store.Directory;
import org.apache.lucene.tests.geo.GeoTestUtil;
import org.apache.lucene.tests.index.RandomIndexWriter;
import org.apache.lucene.tests.search.FixedBitSetCollector;
import org.apache.lucene.tests.util.LuceneTestCase;
import org.apache.lucene.tests.util.TestUtil;
import org.apache.lucene.util.DocIdSetBuilder;
import org.apache.lucene.util.FixedBitSet;
import org.apache.lucene.util.IOUtils;
import org.apache.lucene.util.NumericUtils;

public class TestGeo3DPoint extends LuceneTestCase {

  protected PlanetModel randomPlanetModel() {
    return RandomPicks.randomFrom(
        random(),
        new PlanetModel[] {PlanetModel.WGS84, PlanetModel.CLARKE_1866, PlanetModel.SPHERE});
  }

  private static Codec getCodec() {
    if (Codec.getDefault().getName().equals("Lucene84")) {
      int maxPointsInLeafNode = TestUtil.nextInt(random(), 16, 2048);
      double maxMBSortInHeap = 3.0 + (3 * random().nextDouble());
      if (VERBOSE) {
        System.out.println(
            "TEST: using Lucene60PointsFormat with maxPointsInLeafNode="
                + maxPointsInLeafNode
                + " and maxMBSortInHeap="
                + maxMBSortInHeap);
      }

      return new FilterCodec("Lucene84", Codec.getDefault()) {
        @Override
        public PointsFormat pointsFormat() {
          return new PointsFormat() {
            @Override
            public PointsWriter fieldsWriter(SegmentWriteState writeState) throws IOException {
              return new Lucene90PointsWriter(writeState, maxPointsInLeafNode, maxMBSortInHeap);
            }

            @Override
            public PointsReader fieldsReader(SegmentReadState readState) throws IOException {
              return new Lucene90PointsReader(readState);
            }
          };
        }
      };
    } else {
      return Codec.getDefault();
    }
  }

  public void testBasic() throws Exception {
    Directory dir = getDirectory();
    IndexWriterConfig iwc = newIndexWriterConfig();
    iwc.setCodec(getCodec());
    IndexWriter w = new IndexWriter(dir, iwc);
    Document doc = new Document();
    PlanetModel planetModel = randomPlanetModel();
    doc.add(new Geo3DPoint("field", planetModel, 50.7345267, -97.5303555));
    w.addDocument(doc);
    IndexReader r = DirectoryReader.open(w);
    // We can't wrap with "exotic" readers because the query must see the BKD3DDVFormat:
    IndexSearcher s = newSearcher(r, false);
    assertEquals(
        1,
        s.search(
                Geo3DPoint.newShapeQuery(
                    "field",
                    GeoCircleFactory.makeGeoCircle(
                        planetModel, toRadians(50), toRadians(-97), Math.PI / 180.)),
                1)
            .totalHits
            .value());
    w.close();
    r.close();
    dir.close();
  }

  private static double toRadians(double degrees) {
    return degrees * Geo3DUtil.RADIANS_PER_DEGREE;
  }

  private static class Cell {
    static final IntSupplier nextCellID =
        new IntSupplier() {
          int counter = 0;

          @Override
          public int getAsInt() {
            return counter++;
          }
        };

    final Cell parent;
    final int cellID;
    final int xMinEnc, xMaxEnc;
    final int yMinEnc, yMaxEnc;
    final int zMinEnc, zMaxEnc;
    final int splitCount;
    final PlanetModel planetModel;

    public Cell(
        Cell parent,
        int xMinEnc,
        int xMaxEnc,
        int yMinEnc,
        int yMaxEnc,
        int zMinEnc,
        int zMaxEnc,
        final PlanetModel planetModel,
        int splitCount) {
      this.parent = parent;
      this.xMinEnc = xMinEnc;
      this.xMaxEnc = xMaxEnc;
      this.yMinEnc = yMinEnc;
      this.yMaxEnc = yMaxEnc;
      this.zMinEnc = zMinEnc;
      this.zMaxEnc = zMaxEnc;
      this.cellID = nextCellID.getAsInt();
      this.splitCount = splitCount;
      this.planetModel = planetModel;
    }

    /** Returns true if the quantized point lies within this cell, inclusive on all bounds. */
    public boolean contains(GeoPoint point) {
      int docX = planetModel.encodeValue(point.x);
      int docY = planetModel.encodeValue(point.y);
      int docZ = planetModel.encodeValue(point.z);

      return docX >= xMinEnc
          && docX <= xMaxEnc
          && docY >= yMinEnc
          && docY <= yMaxEnc
          && docZ >= zMinEnc
          && docZ <= zMaxEnc;
    }

    @Override
    public String toString() {
      return "cell="
          + cellID
          + (parent == null ? "" : " parentCellID=" + parent.cellID)
          + " x: "
          + xMinEnc
          + " TO "
          + xMaxEnc
          + ", y: "
          + yMinEnc
          + " TO "
          + yMaxEnc
          + ", z: "
          + zMinEnc
          + " TO "
          + zMaxEnc
          + ", splits: "
          + splitCount;
    }
  }

  private static double quantize(double xyzValue, final PlanetModel planetModel) {
    return planetModel.decodeValue(planetModel.encodeValue(xyzValue));
  }

  private static GeoPoint quantize(GeoPoint point, final PlanetModel planetModel) {
    return new GeoPoint(
        quantize(point.x, planetModel),
        quantize(point.y, planetModel),
        quantize(point.z, planetModel));
  }

  /** Tests consistency of GeoArea.getRelationship vs GeoShape.isWithin */
  public void testGeo3DRelations() throws Exception {

    int numDocs = atLeast(200);
    if (VERBOSE) {
      System.out.println("TEST: " + numDocs + " docs");
    }

    GeoPoint[] docs = new GeoPoint[numDocs];
    GeoPoint[] unquantizedDocs = new GeoPoint[numDocs];
    PlanetModel planetModel = PlanetModel.CLARKE_1866;
    for (int docID = 0; docID < numDocs; docID++) {
      unquantizedDocs[docID] =
          new GeoPoint(
              planetModel,
              toRadians(GeoTestUtil.nextLatitude()),
              toRadians(GeoTestUtil.nextLongitude()));
      docs[docID] = quantize(unquantizedDocs[docID], planetModel);
      if (VERBOSE) {
        System.out.println(
            "  doc=" + docID + ": " + docs[docID] + "; unquantized: " + unquantizedDocs[docID]);
      }
    }

    int iters = atLeast(10);

    int recurseDepth = RandomNumbers.randomIntBetween(random(), 5, 15);

    for (int iter = 0; iter < iters; iter++) {
      GeoShape shape = randomShape(planetModel);

      StringWriter sw = new StringWriter();
      PrintWriter log = new PrintWriter(sw, true);

      if (VERBOSE) {
        log.println("TEST: iter=" + iter + " shape=" + shape);
      }

      XYZBounds bounds = new XYZBounds();
      shape.getBounds(bounds);

      // Start with the root cell that fully contains the shape:
      Cell root =
          new Cell(
              null,
              encodeValueLenient(bounds.getMinimumX(), planetModel),
              encodeValueLenient(bounds.getMaximumX(), planetModel),
              encodeValueLenient(bounds.getMinimumY(), planetModel),
              encodeValueLenient(bounds.getMaximumY(), planetModel),
              encodeValueLenient(bounds.getMinimumZ(), planetModel),
              encodeValueLenient(bounds.getMaximumZ(), planetModel),
              planetModel,
              0);

      if (VERBOSE) {
        log.println("  root cell: " + root);
      }

      // make sure the root cell (XYZBounds) does in fact contain all points that the shape contains
      {
        boolean fail = false;
        for (int docID = 0; docID < numDocs; docID++) {
          if (root.contains(docs[docID]) == false) {
            boolean expected = shape.isWithin(unquantizedDocs[docID]);
            if (expected) {
              log.println(
                  "    doc="
                      + docID
                      + " is contained by shape but is outside the returned XYZBounds");
              log.println("      unquantized=" + unquantizedDocs[docID]);
              log.println("      quantized=" + docs[docID]);
              fail = true;
            }
          }
        }

        if (fail) {
          log.println("  shape=" + shape);
          log.println("  bounds=" + bounds);
          System.out.print(sw.toString());
          fail("invalid bounds for shape=" + shape);
        }
      }

      List<Cell> queue = new ArrayList<>();
      queue.add(root);
      Set<Integer> hits = new HashSet<>();

      while (queue.size() > 0) {
        Cell cell = queue.get(queue.size() - 1);
        queue.remove(queue.size() - 1);
        if (VERBOSE) {
          log.println("  cycle: " + cell + " queue.size()=" + queue.size());
        }

        if (random().nextInt(10) == 7 || cell.splitCount > recurseDepth) {
          if (VERBOSE) {
            log.println("    leaf");
          }
          // Leaf cell: brute force check all docs that fall within this cell:
          for (int docID = 0; docID < numDocs; docID++) {
            GeoPoint point = docs[docID];
            GeoPoint mappedPoint = unquantizedDocs[docID];
            boolean pointWithinShape = shape.isWithin(point);
            boolean mappedPointWithinShape = shape.isWithin(mappedPoint);
            if (cell.contains(point)) {
              if (mappedPointWithinShape) {
                if (VERBOSE) {
                  log.println(
                      "    check doc="
                          + docID
                          + ": match!  Actual quantized point within: "
                          + pointWithinShape);
                }
                hits.add(docID);
              } else {
                if (VERBOSE) {
                  log.println(
                      "    check doc="
                          + docID
                          + ": no match.  Quantized point within: "
                          + pointWithinShape);
                }
              }
            }
          }
        } else {
          GeoArea xyzSolid =
              GeoAreaFactory.makeGeoArea(
                  planetModel,
                  Geo3DUtil.decodeValueFloor(cell.xMinEnc, planetModel),
                  Geo3DUtil.decodeValueCeil(cell.xMaxEnc, planetModel),
                  Geo3DUtil.decodeValueFloor(cell.yMinEnc, planetModel),
                  Geo3DUtil.decodeValueCeil(cell.yMaxEnc, planetModel),
                  Geo3DUtil.decodeValueFloor(cell.zMinEnc, planetModel),
                  Geo3DUtil.decodeValueCeil(cell.zMaxEnc, planetModel));

          if (VERBOSE) {
            log.println(
                "    minx="
                    + Geo3DUtil.decodeValueFloor(cell.xMinEnc, planetModel)
                    + " maxx="
                    + Geo3DUtil.decodeValueCeil(cell.xMaxEnc, planetModel)
                    + " miny="
                    + Geo3DUtil.decodeValueFloor(cell.yMinEnc, planetModel)
                    + " maxy="
                    + Geo3DUtil.decodeValueCeil(cell.yMaxEnc, planetModel)
                    + " minz="
                    + Geo3DUtil.decodeValueFloor(cell.zMinEnc, planetModel)
                    + " maxz="
                    + Geo3DUtil.decodeValueCeil(cell.zMaxEnc, planetModel));
          }

          switch (xyzSolid.getRelationship(shape)) {
            case GeoArea.CONTAINS:
              // Shape fully contains the cell: blindly add all docs in this cell:
              if (VERBOSE) {
                log.println("    GeoArea.CONTAINS: now addAll");
              }
              for (int docID = 0; docID < numDocs; docID++) {
                if (cell.contains(docs[docID])) {
                  if (VERBOSE) {
                    log.println("    addAll doc=" + docID);
                  }
                  hits.add(docID);
                }
              }
              continue;
            case GeoArea.OVERLAPS:
              if (VERBOSE) {
                log.println("    GeoArea.OVERLAPS: keep splitting");
              }
              // They do overlap but neither contains the other:
              // log.println("    crosses1");
              break;
            case GeoArea.WITHIN:
              if (VERBOSE) {
                log.println("    GeoArea.WITHIN: keep splitting");
              }
              // Cell fully contains the shape:
              // log.println("    crosses2");
              break;
            case GeoArea.DISJOINT:
              // They do not overlap at all: don't recurse on this cell
              // log.println("    outside");
              if (VERBOSE) {
                log.println("    GeoArea.DISJOINT: drop this cell");
                for (int docID = 0; docID < numDocs; docID++) {
                  if (cell.contains(docs[docID])) {
                    log.println("    skip doc=" + docID);
                  }
                }
              }
              continue;
            default:
              assert false;
          }

          // Randomly split:
          switch (random().nextInt(3)) {
            case 0:
              // Split on X:
              {
                int splitValue =
                    RandomNumbers.randomIntBetween(random(), cell.xMinEnc, cell.xMaxEnc);
                if (VERBOSE) {
                  log.println("    now split on x=" + splitValue);
                }
                Cell cell1 =
                    new Cell(
                        cell,
                        cell.xMinEnc,
                        splitValue,
                        cell.yMinEnc,
                        cell.yMaxEnc,
                        cell.zMinEnc,
                        cell.zMaxEnc,
                        planetModel,
                        cell.splitCount + 1);
                Cell cell2 =
                    new Cell(
                        cell,
                        splitValue,
                        cell.xMaxEnc,
                        cell.yMinEnc,
                        cell.yMaxEnc,
                        cell.zMinEnc,
                        cell.zMaxEnc,
                        planetModel,
                        cell.splitCount + 1);
                if (VERBOSE) {
                  log.println("    split cell1: " + cell1);
                  log.println("    split cell2: " + cell2);
                }
                queue.add(cell1);
                queue.add(cell2);
              }
              break;

            case 1:
              // Split on Y:
              {
                int splitValue =
                    RandomNumbers.randomIntBetween(random(), cell.yMinEnc, cell.yMaxEnc);
                if (VERBOSE) {
                  log.println("    now split on y=" + splitValue);
                }
                Cell cell1 =
                    new Cell(
                        cell,
                        cell.xMinEnc,
                        cell.xMaxEnc,
                        cell.yMinEnc,
                        splitValue,
                        cell.zMinEnc,
                        cell.zMaxEnc,
                        planetModel,
                        cell.splitCount + 1);
                Cell cell2 =
                    new Cell(
                        cell,
                        cell.xMinEnc,
                        cell.xMaxEnc,
                        splitValue,
                        cell.yMaxEnc,
                        cell.zMinEnc,
                        cell.zMaxEnc,
                        planetModel,
                        cell.splitCount + 1);
                if (VERBOSE) {
                  log.println("    split cell1: " + cell1);
                  log.println("    split cell2: " + cell2);
                }
                queue.add(cell1);
                queue.add(cell2);
              }
              break;

            case 2:
              // Split on Z:
              {
                int splitValue =
                    RandomNumbers.randomIntBetween(random(), cell.zMinEnc, cell.zMaxEnc);
                if (VERBOSE) {
                  log.println("    now split on z=" + splitValue);
                }
                Cell cell1 =
                    new Cell(
                        cell,
                        cell.xMinEnc,
                        cell.xMaxEnc,
                        cell.yMinEnc,
                        cell.yMaxEnc,
                        cell.zMinEnc,
                        splitValue,
                        planetModel,
                        cell.splitCount + 1);
                Cell cell2 =
                    new Cell(
                        cell,
                        cell.xMinEnc,
                        cell.xMaxEnc,
                        cell.yMinEnc,
                        cell.yMaxEnc,
                        splitValue,
                        cell.zMaxEnc,
                        planetModel,
                        cell.splitCount + 1);
                if (VERBOSE) {
                  log.println("    split cell1: " + cell1);
                  log.println("    split cell2: " + cell2);
                }
                queue.add(cell1);
                queue.add(cell2);
              }
              break;
          }
        }
      }

      if (VERBOSE) {
        log.println("  " + hits.size() + " hits");
      }

      // Done matching, now verify:
      boolean fail = false;
      for (int docID = 0; docID < numDocs; docID++) {
        GeoPoint point = docs[docID];
        GeoPoint mappedPoint = unquantizedDocs[docID];
        boolean expected = shape.isWithin(mappedPoint);
        boolean actual = hits.contains(docID);
        if (actual != expected) {
          if (actual) {
            log.println("doc=" + docID + " should not have matched but did");
          } else {
            log.println("doc=" + docID + " should match but did not");
          }
          log.println("  point=" + point);
          log.println("  mappedPoint=" + mappedPoint);
          fail = true;
        }
      }

      if (fail) {
        System.out.print(sw.toString());
        fail("invalid hits for shape=" + shape);
      }
    }
  }

  public void testRandomTiny() throws Exception {
    // Make sure single-leaf-node case is OK:
    doTestRandom(10);
  }

  // TODO: incredibly slow
  @Nightly
  public void testRandomMedium() throws Exception {
    doTestRandom(1000);
  }

  @Nightly
  public void testRandomBig() throws Exception {
    doTestRandom(50000);
  }

  private void doTestRandom(int count) throws Exception {
    int numPoints = atLeast(count);

    if (VERBOSE) {
      System.err.println("TEST: numPoints=" + numPoints);
    }

    double[] lats = new double[numPoints];
    double[] lons = new double[numPoints];

    boolean haveRealDoc = false;

    Random random = random();
    for (int docID = 0; docID < numPoints; docID++) {
      int x = random.nextInt(20);
      if (x == 17) {
        // Some docs don't have a point:
        lats[docID] = Double.NaN;
        if (VERBOSE) {
          System.err.println("  doc=" + docID + " is missing");
        }
        continue;
      }

      if (docID > 0 && x < 3 && haveRealDoc) {
        int oldDocID;
        while (true) {
          oldDocID = random.nextInt(docID);
          if (Double.isNaN(lats[oldDocID]) == false) {
            break;
          }
        }

        if (x == 0) {
          // Identical lat to old point
          lats[docID] = lats[oldDocID];
          lons[docID] = GeoTestUtil.nextLongitude();
          if (VERBOSE) {
            System.err.println(
                "  doc="
                    + docID
                    + " lat="
                    + lats[docID]
                    + " lon="
                    + lons[docID]
                    + " (same lat as doc="
                    + oldDocID
                    + ")");
          }
        } else if (x == 1) {
          // Identical lon to old point
          lats[docID] = GeoTestUtil.nextLatitude();
          lons[docID] = lons[oldDocID];
          if (VERBOSE) {
            System.err.println(
                "  doc="
                    + docID
                    + " lat="
                    + lats[docID]
                    + " lon="
                    + lons[docID]
                    + " (same lon as doc="
                    + oldDocID
                    + ")");
          }
        } else {
          assert x == 2;
          // Fully identical point:
          lats[docID] = lats[oldDocID];
          lons[docID] = lons[oldDocID];
          if (VERBOSE) {
            System.err.println(
                "  doc="
                    + docID
                    + " lat="
                    + lats[docID]
                    + " lon="
                    + lons[docID]
                    + " (same lat/lon as doc="
                    + oldDocID
                    + ")");
          }
        }
      } else {
        lats[docID] = GeoTestUtil.nextLatitude();
        lons[docID] = GeoTestUtil.nextLongitude();
        haveRealDoc = true;
        if (VERBOSE) {
          System.err.println("  doc=" + docID + " lat=" + lats[docID] + " lon=" + lons[docID]);
        }
      }
    }

    verify(lats, lons, randomPlanetModel());
  }

  public void testPolygonOrdering() {
    final double[] lats =
        new double[] {
          51.204382859999996,
          50.89947531437482,
          50.8093624806861,
          50.8093624806861,
          50.89947531437482,
          51.204382859999996,
          51.51015366140113,
          51.59953838204167,
          51.59953838204167,
          51.51015366140113,
          51.204382859999996
        };
    final double[] lons =
        new double[] {
          0.8747711978759765,
          0.6509219832137298,
          0.35960265165247807,
          0.10290284834752167,
          -0.18841648321373008,
          -0.41226569787597667,
          -0.18960465285650027,
          0.10285893781346236,
          0.35964656218653757,
          0.6521101528565002,
          0.8747711978759765
        };
    final Query q =
        Geo3DPoint.newPolygonQuery("point", randomPlanetModel(), new Polygon(lats, lons));
    // System.out.println(q);
    assertTrue(!q.toString().contains("GeoConcavePolygon"));
  }

  private static Query random3DQuery(final String field, final PlanetModel planetModel) {
    while (true) {
      final int shapeType = random().nextInt(5);
      switch (shapeType) {
        case 4:
          {
            // Large polygons
            final boolean isClockwise = random().nextDouble() < 0.5;
            try {
              final Query q =
                  Geo3DPoint.newLargePolygonQuery(
                      field,
                      planetModel,
                      makePoly(
                          planetModel,
                          new GeoPoint(
                              planetModel,
                              toRadians(GeoTestUtil.nextLatitude()),
                              toRadians(GeoTestUtil.nextLongitude())),
                          isClockwise,
                          true));
              // System.err.println("Generated: "+q);
              // assertTrue(false);
              return q;
            } catch (IllegalArgumentException _) {
              continue;
            }
          }

        case 0:
          {
            // Polygons
            final boolean isClockwise = random().nextDouble() < 0.5;
            try {
              final Query q =
                  Geo3DPoint.newPolygonQuery(
                      field,
                      planetModel,
                      makePoly(
                          planetModel,
                          new GeoPoint(
                              planetModel,
                              toRadians(GeoTestUtil.nextLatitude()),
                              toRadians(GeoTestUtil.nextLongitude())),
                          isClockwise,
                          true));
              // System.err.println("Generated: "+q);
              // assertTrue(false);
              return q;
            } catch (IllegalArgumentException _) {
              continue;
            }
          }

        case 1:
          {
            // Circles
            final double widthMeters =
                random().nextDouble() * Math.PI * planetModel.getMeanRadius();
            try {
              return Geo3DPoint.newDistanceQuery(
                  field,
                  planetModel,
                  GeoTestUtil.nextLatitude(),
                  GeoTestUtil.nextLongitude(),
                  widthMeters);
            } catch (IllegalArgumentException _) {
              continue;
            }
          }

        case 2:
          {
            // Rectangles
            final Rectangle r = GeoTestUtil.nextBox();
            try {
              return Geo3DPoint.newBoxQuery(
                  field, planetModel, r.minLat, r.maxLat, r.minLon, r.maxLon);
            } catch (IllegalArgumentException _) {
              continue;
            }
          }

        case 3:
          {
            // Paths
            // TBD: Need to rework generation to be realistic
            final int pointCount = random().nextInt(5) + 1;
            final double width =
                random().nextDouble() * Math.PI * 0.5 * planetModel.getMeanRadius();
            final double[] latitudes = new double[pointCount];
            final double[] longitudes = new double[pointCount];
            for (int i = 0; i < pointCount; i++) {
              latitudes[i] = GeoTestUtil.nextLatitude();
              longitudes[i] = GeoTestUtil.nextLongitude();
            }
            try {
              return Geo3DPoint.newPathQuery(field, latitudes, longitudes, width, planetModel);
            } catch (IllegalArgumentException _) {
              // This is what happens when we create a shape that is invalid.  Although it is
              // conceivable that there are cases where
              // the exception is thrown incorrectly, we aren't going to be able to do that in this
              // random test.
              continue;
            }
          }

        default:
          throw new IllegalStateException("Unexpected shape type");
      }
    }
  }

  // Poached from Geo3dRptTest.randomShape:
  private static GeoShape randomShape(final PlanetModel planetModel) {
    while (true) {
      final int shapeType = random().nextInt(4);
      switch (shapeType) {
        case 0:
          {
            // Polygons
            final int vertexCount = random().nextInt(3) + 3;
            final List<GeoPoint> geoPoints = new ArrayList<>();
            while (geoPoints.size() < vertexCount) {
              final GeoPoint gPt =
                  new GeoPoint(
                      planetModel,
                      toRadians(GeoTestUtil.nextLatitude()),
                      toRadians(GeoTestUtil.nextLongitude()));
              geoPoints.add(gPt);
            }
            try {
              final GeoShape rval = GeoPolygonFactory.makeGeoPolygon(planetModel, geoPoints);
              if (rval == null) {
                // Degenerate polygon
                continue;
              }
              return rval;
            } catch (IllegalArgumentException _) {
              // This is what happens when we create a shape that is invalid.  Although it is
              // conceivable that there are cases where
              // the exception is thrown incorrectly, we aren't going to be able to do that in this
              // random test.
              continue;
            }
          }

        case 1:
          {
            // Circles

            double lat = toRadians(GeoTestUtil.nextLatitude());
            double lon = toRadians(GeoTestUtil.nextLongitude());

            double angle = random().nextDouble() * Math.PI / 2.0;

            try {
              return GeoCircleFactory.makeGeoCircle(planetModel, lat, lon, angle);
            } catch (IllegalArgumentException _) {
              // angle is too small; try again:
              continue;
            }
          }

        case 2:
          {
            // Rectangles
            double lat0 = toRadians(GeoTestUtil.nextLatitude());
            double lat1 = toRadians(GeoTestUtil.nextLatitude());
            if (lat1 < lat0) {
              double x = lat0;
              lat0 = lat1;
              lat1 = x;
            }
            double lon0 = toRadians(GeoTestUtil.nextLongitude());
            double lon1 = toRadians(GeoTestUtil.nextLongitude());
            if (lon1 < lon0) {
              double x = lon0;
              lon0 = lon1;
              lon1 = x;
            }

            return GeoBBoxFactory.makeGeoBBox(planetModel, lat1, lat0, lon0, lon1);
          }

        case 3:
          {
            // Paths
            final int pointCount = random().nextInt(5) + 1;
            final double width = toRadians(random().nextInt(89) + 1);
            final GeoPoint[] points = new GeoPoint[pointCount];
            for (int i = 0; i < pointCount; i++) {
              points[i] =
                  new GeoPoint(
                      planetModel,
                      toRadians(GeoTestUtil.nextLatitude()),
                      toRadians(GeoTestUtil.nextLongitude()));
            }
            try {
              return GeoPathFactory.makeGeoPath(planetModel, width, points);
            } catch (IllegalArgumentException _) {
              // This is what happens when we create a shape that is invalid.  Although it is
              // conceivable that there are cases where
              // the exception is thrown incorrectly, we aren't going to be able to do that in this
              // random test.
              continue;
            }
          }

        default:
          throw new IllegalStateException("Unexpected shape type");
      }
    }
  }

  private static void verify(double[] lats, double[] lons, final PlanetModel planetModel)
      throws Exception {
    IndexWriterConfig iwc = newIndexWriterConfig();

    GeoPoint[] points = new GeoPoint[lats.length];
    GeoPoint[] unquantizedPoints = new GeoPoint[lats.length];

    // Pre-quantize all lat/lons:
    for (int i = 0; i < lats.length; i++) {
      if (Double.isNaN(lats[i]) == false) {
        // System.out.println("lats[" + i + "] = " + lats[i]);
        unquantizedPoints[i] = new GeoPoint(planetModel, toRadians(lats[i]), toRadians(lons[i]));
        points[i] = quantize(unquantizedPoints[i], planetModel);
      }
    }

    // Else we can get O(N^2) merging:
    int mbd = iwc.getMaxBufferedDocs();
    if (mbd != -1 && mbd < points.length / 100) {
      iwc.setMaxBufferedDocs(points.length / 100);
    }
    iwc.setCodec(getCodec());
    Directory dir;
    if (points.length > 100000) {
      dir = newFSDirectory(createTempDir("TestBKDTree"));
    } else {
      dir = getDirectory();
    }
    Set<Integer> deleted = new HashSet<>();
    // RandomIndexWriter is too slow here:
    IndexWriter w = new IndexWriter(dir, iwc);
    for (int id = 0; id < points.length; id++) {
      Document doc = new Document();
      doc.add(new StringField("id", "" + id, Field.Store.NO));
      doc.add(new NumericDocValuesField("id", id));
      GeoPoint point = points[id];
      if (point != null) {
        doc.add(new Geo3DPoint("point", planetModel, point.x, point.y, point.z));
      }
      w.addDocument(doc);
      if (id > 0 && random().nextInt(100) == 42) {
        int idToDelete = random().nextInt(id);
        w.deleteDocuments(new Term("id", "" + idToDelete));
        deleted.add(idToDelete);
        if (VERBOSE) {
          System.err.println("  delete id=" + idToDelete);
        }
      }
    }
    if (random().nextBoolean()) {
      w.forceMerge(1);
    }
    final IndexReader r = DirectoryReader.open(w);
    if (VERBOSE) {
      System.out.println("TEST: using reader " + r);
    }
    w.close();

    // We can't wrap with "exotic" readers because the geo3d query must see the Geo3DDVFormat:
    IndexSearcher s = newSearcher(r, false);

    final int iters = atLeast(100);

    for (int iter = 0; iter < iters; iter++) {

      /*
      GeoShape shape = randomShape();

      if (VERBOSE) {
        System.err.println("\nTEST: iter=" + iter + " shape="+shape);
      }
      */

      Query query =
          random3DQuery("point", planetModel); // Geo3DPoint.newShapeQuery("point", shape);

      if (VERBOSE) {
        System.err.println("  using query: " + query);
      }

      final FixedBitSet hits = s.search(query, FixedBitSetCollector.createManager(r.maxDoc()));

      if (VERBOSE) {
        System.err.println("  hitCount: " + hits.cardinality());
      }

      NumericDocValues docIDToID = MultiDocValues.getNumericValues(r, "id");

      for (int docID = 0; docID < r.maxDoc(); docID++) {
        assertEquals(docID, docIDToID.nextDoc());
        int id = (int) docIDToID.longValue();
        GeoPoint point = points[id];
        GeoPoint unquantizedPoint = unquantizedPoints[id];
        if (point != null && unquantizedPoint != null) {
          GeoShape shape = ((PointInGeo3DShapeQuery) query).getShape();
          XYZBounds bounds = new XYZBounds();
          shape.getBounds(bounds);
          XYZSolid solid =
              XYZSolidFactory.makeXYZSolid(
                  planetModel,
                  bounds.getMinimumX(),
                  bounds.getMaximumX(),
                  bounds.getMinimumY(),
                  bounds.getMaximumY(),
                  bounds.getMinimumZ(),
                  bounds.getMaximumZ());

          boolean expected = ((deleted.contains(id) == false) && shape.isWithin(point));
          if (hits.get(docID) != expected) {
            StringBuilder b = new StringBuilder();
            if (expected) {
              b.append("FAIL: id=" + id + " should have matched but did not\n");
            } else {
              b.append("FAIL: id=" + id + " should not have matched but did\n");
            }
            b.append("  shape=" + shape + "\n");
            b.append("  bounds=" + bounds + "\n");
            b.append(
                "  world bounds=("
                    + " minX="
                    + planetModel.getMinimumXValue()
                    + " maxX="
                    + planetModel.getMaximumXValue()
                    + " minY="
                    + planetModel.getMinimumYValue()
                    + " maxY="
                    + planetModel.getMaximumYValue()
                    + " minZ="
                    + planetModel.getMinimumZValue()
                    + " maxZ="
                    + planetModel.getMaximumZValue()
                    + "\n");
            b.append(
                "  quantized point="
                    + point
                    + " within shape? "
                    + shape.isWithin(point)
                    + " within bounds? "
                    + solid.isWithin(point)
                    + "\n");
            b.append(
                "  unquantized point="
                    + unquantizedPoint
                    + " within shape? "
                    + shape.isWithin(unquantizedPoint)
                    + " within bounds? "
                    + solid.isWithin(unquantizedPoint)
                    + "\n");
            b.append("  docID=" + docID + " deleted?=" + deleted.contains(id) + "\n");
            b.append("  query=" + query + "\n");
            b.append(
                "  explanation:\n    "
                    + explain("point", shape, point, unquantizedPoint, r, docID)
                        .replace("\n", "\n  "));
            fail(b.toString());
          }
        } else {
          assertFalse(hits.get(docID));
        }
      }
    }

    IOUtils.close(r, dir);
  }

  public void testToString() {
    // Don't compare entire strings because Java 9 and Java 8 have slightly different values
    Geo3DPoint point = new Geo3DPoint("point", randomPlanetModel(), 44.244272, 7.769736);
    final String stringToCompare = "Geo3DPoint <point: x=";
    assertEquals(stringToCompare, point.toString().substring(0, stringToCompare.length()));
  }

  public void testShapeQueryToString() {
    // Don't compare entire strings because Java 9 and Java 8 have slightly different values
    final String stringToCompare =
        "PointInGeo3DShapeQuery: field=point: Shape: GeoStandardCircle: {planetmodel=PlanetModel.WGS84, center=[lat=0.7";
    assertEquals(
        stringToCompare,
        Geo3DPoint.newShapeQuery(
                "point",
                GeoCircleFactory.makeGeoCircle(
                    PlanetModel.WGS84, toRadians(44.244272), toRadians(7.769736), 0.1))
            .toString()
            .substring(0, stringToCompare.length()));
  }

  private static Directory getDirectory() {
    return newDirectory();
  }

  public void testEquals() {
    PlanetModel planetModel = randomPlanetModel();
    GeoShape shape = randomShape(planetModel);
    Query q = Geo3DPoint.newShapeQuery("point", shape);
    assertEquals(q, Geo3DPoint.newShapeQuery("point", shape));
    assertFalse(q.equals(Geo3DPoint.newShapeQuery("point2", shape)));

    // make a different random shape:
    GeoShape shape2;
    do {
      shape2 = randomShape(planetModel);
    } while (shape.equals(shape2));

    assertFalse(q.equals(Geo3DPoint.newShapeQuery("point", shape2)));
  }

  public void testComplexPolygons() {
    final PlanetModel pm = randomPlanetModel();
    // Pick a random pole
    final GeoPoint randomPole =
        new GeoPoint(
            pm, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
    int iters = atLeast(100);
    for (int i = 0; i < iters; i++) {
      // Create a polygon that's less than 180 degrees
      makePoly(pm, randomPole, true, true);
    }
    iters = atLeast(100);
    for (int i = 0; i < iters; i++) {
      // Create a polygon that's greater than 180 degrees
      makePoly(pm, randomPole, false, true);
    }
  }

  private static final double MINIMUM_EDGE_ANGLE = toRadians(5.0);
  private static final double MINIMUM_ARC_ANGLE = toRadians(1.0);

  /**
   * Cook up a random Polygon that makes sense, with possible nested polygon within. This is part of
   * testing more complex polygons with nested holes. Picking random points doesn't do it because
   * it's almost impossible to come up with nested ones of the proper clockwise/counterclockwise
   * rotation that way.
   */
  protected static Polygon makePoly(
      final PlanetModel pm,
      final GeoPoint pole,
      final boolean clockwiseDesired,
      final boolean createHoles) {

    // Polygon edges will be arranged around the provided pole, and holes will each have a pole
    // selected within the parent
    // polygon.
    final int pointCount = TestUtil.nextInt(random(), 3, 10);
    // The point angles we pick next.  The only requirement is that they are not all on one side of
    // the pole.
    // We arrange that by picking the next point within what's left of the remaining angle, but
    // never more than 180 degrees,
    // and never less than what is needed to insure that the remaining point choices are less than
    // 180 degrees always.
    // These are all picked in the context of the pole,
    final double[] angles = new double[pointCount];
    final double[] arcDistance = new double[pointCount];
    // Pick a set of points
    while (true) {
      double accumulatedAngle = 0.0;
      for (int i = 0; i < pointCount; i++) {
        final int remainingEdgeCount = pointCount - i;
        final double remainingAngle = 2.0 * Math.PI - accumulatedAngle;
        if (remainingEdgeCount == 1) {
          angles[i] = remainingAngle;
        } else {
          // The maximum angle is 180 degrees, or what's left when you give a minimal amount to each
          // edge.
          double maximumAngle = remainingAngle - (remainingEdgeCount - 1) * MINIMUM_EDGE_ANGLE;
          if (maximumAngle > Math.PI) {
            maximumAngle = Math.PI;
          }
          // The minimum angle is MINIMUM_EDGE_ANGLE, or enough to be sure nobody afterwards needs
          // more than 180 degrees.  And since we have three points to start with, we already know
          // that.
          final double minimumAngle = MINIMUM_EDGE_ANGLE;
          // Pick the angle
          final double angle = random().nextDouble() * (maximumAngle - minimumAngle) + minimumAngle;
          angles[i] = angle;
          accumulatedAngle += angle;
        }
        // Pick the arc distance randomly; not quite the full range though
        arcDistance[i] =
            random().nextDouble() * (Math.PI * 0.5 - MINIMUM_ARC_ANGLE) + MINIMUM_ARC_ANGLE;
      }
      if (clockwiseDesired) {
        // Reverse the signs
        for (int i = 0; i < pointCount; i++) {
          angles[i] = -angles[i];
        }
      }

      // Now, use the pole's information plus angles and arcs to create GeoPoints in the right
      // order.
      final List<GeoPoint> polyPoints = convertToPoints(pm, pole, angles, arcDistance);

      // Next, do some holes.  No more than 2 of these.  The poles for holes must always be within
      // the polygon, so we're going to use Geo3D to help us select those given the points we just
      // made.

      final List<Polygon> holeList = new ArrayList<>();

      /* Hole logic is broken and needs rethinking

      final int holeCount = createHoles ? TestUtil.nextInt(random(), 0, 2) : 0;

      // Create the geo3d polygon, so we can test out our poles.
      final GeoPolygon poly;
      try {
        poly = GeoPolygonFactory.makeGeoPolygon(pm, polyPoints, null);
      } catch (IllegalArgumentException e) {
        // This is what happens when three adjacent points are colinear, so try again.
        continue;
      }

      for (int i = 0; i < holeCount; i++) {
        // Choose a pole.  The poly has to be within the polygon, but it also cannot be on the polygon edge.
        // If we can't find a good pole we have to give it up and not do the hole.
        for (int k = 0; k < 500; k++) {
          final GeoPoint poleChoice = new GeoPoint(pm, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
          if (!poly.isWithin(poleChoice)) {
            continue;
          }
          // We have a pole within the polygon.  Now try 100 times to build a polygon that does not intersect the outside ring.
          // After that we give up and pick a new pole.
          boolean foundOne = false;
          for (int j = 0; j < 100; j++) {
            final Polygon insidePoly = makePoly(pm, poleChoice, !clockwiseDesired, false);
            // Verify that the inside polygon is OK.  If not, discard and repeat.
            if (!verifyPolygon(pm, insidePoly, poly)) {
              continue;
            }
            holeList.add(insidePoly);
            foundOne = true;
          }
          if (foundOne) {
            break;
          }
        }
      }
      */

      final Polygon[] holes = holeList.toArray(new Polygon[0]);

      // Finally, build the polygon and return it
      final double[] lats = new double[polyPoints.size() + 1];
      final double[] lons = new double[polyPoints.size() + 1];

      for (int i = 0; i < polyPoints.size(); i++) {
        lats[i] = polyPoints.get(i).getLatitude() * 180.0 / Math.PI;
        lons[i] = polyPoints.get(i).getLongitude() * 180.0 / Math.PI;
      }
      lats[polyPoints.size()] = lats[0];
      lons[polyPoints.size()] = lons[0];
      return new Polygon(lats, lons, holes);
    }
  }

  protected static List<GeoPoint> convertToPoints(
      final PlanetModel pm,
      final GeoPoint pole,
      final double[] angles,
      final double[] arcDistances) {
    // To do the point rotations, we need the sine and cosine of the pole latitude and longitude.
    // Get it here for performance.
    final double sinLatitude = Math.sin(pole.getLatitude());
    final double cosLatitude = Math.cos(pole.getLatitude());
    final double sinLongitude = Math.sin(pole.getLongitude());
    final double cosLongitude = Math.cos(pole.getLongitude());
    final List<GeoPoint> rval = new ArrayList<>();
    for (int i = 0; i < angles.length; i++) {
      rval.add(
          createPoint(
              pm,
              angles[i],
              arcDistances[i],
              sinLatitude,
              cosLatitude,
              sinLongitude,
              cosLongitude));
    }
    return rval;
  }

  protected static GeoPoint createPoint(
      final PlanetModel pm,
      final double angle,
      final double arcDistance,
      final double sinLatitude,
      final double cosLatitude,
      final double sinLongitude,
      final double cosLongitude) {
    // From the angle and arc distance, convert to (x,y,z) in unit space.
    // We want the perspective to be looking down the x axis.  The "angle" measurement is thus in
    // the Y-Z plane. The arcdistance is in X.
    final double x = Math.cos(arcDistance);
    final double yzScale = Math.sin(arcDistance);
    final double y = Math.cos(angle) * yzScale;
    final double z = Math.sin(angle) * yzScale;
    // Now, rotate coordinates so that we shift everything from pole = x-axis to actual coordinates.
    // This transformation should take the point (1,0,0) and transform it to the pole's actual
    // (x,y,z) coordinates.
    // Coordinate rotation formula:
    // x1 = x0 cos T - y0 sin T
    // y1 = x0 sin T + y0 cos T
    // We're in essence undoing the following transformation (from GeoPolygonFactory):
    // x1 = x0 cos az + y0 sin az
    // y1 = - x0 sin az + y0 cos az
    // z1 = z0
    // x2 = x1 cos al + z1 sin al
    // y2 = y1
    // z2 = - x1 sin al + z1 cos al
    // So, we reverse the order of the transformations, AND we transform backwards.
    // Transforming backwards means using these identities: sin(-angle) = -sin(angle), cos(-angle) =
    // cos(angle)
    // So:
    // x1 = x0 cos al - z0 sin al
    // y1 = y0
    // z1 = x0 sin al + z0 cos al
    // x2 = x1 cos az - y1 sin az
    // y2 = x1 sin az + y1 cos az
    // z2 = z1
    final double x1 = x * cosLatitude - z * sinLatitude;
    final double y1 = y;
    final double z1 = x * sinLatitude + z * cosLatitude;
    final double x2 = x1 * cosLongitude - y1 * sinLongitude;
    final double y2 = x1 * sinLongitude + y1 * cosLongitude;
    final double z2 = z1;
    // Scale to put the point on the surface
    return pm.createSurfacePoint(x2, y2, z2);
  }

  protected static boolean verifyPolygon(
      final PlanetModel pm, final Polygon polygon, final GeoPolygon outsidePolygon) {
    // Each point in the new poly should be inside the outside poly, and each edge should not
    // intersect the outside poly edge
    final double[] lats = polygon.getPolyLats();
    final double[] lons = polygon.getPolyLons();
    final List<GeoPoint> polyPoints = new ArrayList<>(lats.length - 1);
    for (int i = 0; i < lats.length - 1; i++) {
      final GeoPoint newPoint = new GeoPoint(pm, toRadians(lats[i]), toRadians(lons[i]));
      if (!outsidePolygon.isWithin(newPoint)) {
        return false;
      }
      polyPoints.add(newPoint);
    }
    // We don't need to construct the world to find intersections -- just the bordering planes.
    for (int planeIndex = 0; planeIndex < polyPoints.size(); planeIndex++) {
      final GeoPoint startPoint = polyPoints.get(planeIndex);
      final GeoPoint endPoint = polyPoints.get(legalIndex(planeIndex + 1, polyPoints.size()));
      final GeoPoint beforeStartPoint =
          polyPoints.get(legalIndex(planeIndex - 1, polyPoints.size()));
      final GeoPoint afterEndPoint = polyPoints.get(legalIndex(planeIndex + 2, polyPoints.size()));
      final SidedPlane beforePlane = new SidedPlane(endPoint, beforeStartPoint, startPoint);
      final SidedPlane afterPlane = new SidedPlane(startPoint, endPoint, afterEndPoint);
      final Plane plane = new Plane(startPoint, endPoint);

      // Check for intersections!!
      if (outsidePolygon.intersects(plane, null, beforePlane, afterPlane)) {
        return false;
      }
    }
    return true;
  }

  protected static int legalIndex(int index, int size) {
    if (index >= size) {
      index -= size;
    }
    if (index < 0) {
      index += size;
    }
    return index;
  }

  public void testEncodeDecodeCeil() throws Exception {
    PlanetModel planetModel = randomPlanetModel();
    // just for testing quantization error
    final double ENCODING_TOLERANCE = planetModel.DECODE;

    int iters = atLeast(10000);
    for (int iter = 0; iter < iters; iter++) {
      GeoPoint point =
          new GeoPoint(
              planetModel,
              toRadians(GeoTestUtil.nextLatitude()),
              toRadians(GeoTestUtil.nextLongitude()));
      double xEnc = planetModel.decodeValue(planetModel.encodeValue(point.x));
      assertEquals(
          "x=" + point.x + " xEnc=" + xEnc + " diff=" + (point.x - xEnc),
          point.x,
          xEnc,
          ENCODING_TOLERANCE);

      double yEnc = planetModel.decodeValue(planetModel.encodeValue(point.y));
      assertEquals(
          "y=" + point.y + " yEnc=" + yEnc + " diff=" + (point.y - yEnc),
          point.y,
          yEnc,
          ENCODING_TOLERANCE);

      double zEnc = planetModel.decodeValue(planetModel.encodeValue(point.z));
      assertEquals(
          "z=" + point.z + " zEnc=" + zEnc + " diff=" + (point.z - zEnc),
          point.z,
          zEnc,
          ENCODING_TOLERANCE);
    }

    // check edge/interesting cases explicitly
    double planetMax = planetModel.getMaximumMagnitude();
    for (double value : new double[] {0.0, -planetMax, planetMax}) {
      assertEquals(
          value, planetModel.decodeValue(planetModel.encodeValue(value)), ENCODING_TOLERANCE);
    }
  }

  /** make sure values always go down: this is important for edge case consistency */
  public void testEncodeDecodeRoundsDown() throws Exception {
    PlanetModel planetModel = randomPlanetModel();
    int iters = atLeast(1000);
    for (int iter = 0; iter < iters; iter++) {
      final double latBase = GeoTestUtil.nextLatitude();
      final double lonBase = GeoTestUtil.nextLongitude();

      // test above the value
      double lat = latBase;
      double lon = lonBase;
      for (int i = 0; i < 1000; i++) {
        lat = Math.min(90, Math.nextUp(lat));
        lon = Math.min(180, Math.nextUp(lon));
        GeoPoint point = new GeoPoint(planetModel, toRadians(lat), toRadians(lon));
        GeoPoint pointEnc =
            new GeoPoint(
                Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.x), planetModel),
                Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.y), planetModel),
                Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.z), planetModel));
        assertTrue(pointEnc.x <= point.x);
        assertTrue(pointEnc.y <= point.y);
        assertTrue(pointEnc.z <= point.z);
      }

      // test below the value
      lat = latBase;
      lon = lonBase;
      for (int i = 0; i < 1000; i++) {
        lat = Math.max(-90, Math.nextDown(lat));
        lon = Math.max(-180, Math.nextDown(lon));
        GeoPoint point = new GeoPoint(planetModel, toRadians(lat), toRadians(lon));
        GeoPoint pointEnc =
            new GeoPoint(
                Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.x), planetModel),
                Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.y), planetModel),
                Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.z), planetModel));
        assertTrue(pointEnc.x <= point.x);
        assertTrue(pointEnc.y <= point.y);
        assertTrue(pointEnc.z <= point.z);
      }
    }
  }

  public void testMinValueQuantization() {
    PlanetModel planetModel = randomPlanetModel();
    int encoded = planetModel.MIN_ENCODED_VALUE;
    double minValue = -planetModel.getMaximumMagnitude();
    // Normal encoding
    double decoded = planetModel.decodeValue(encoded);
    assertEquals(minValue, decoded, 0d);
    assertEquals(encoded, planetModel.encodeValue(decoded));
    // Encoding floor
    double decodedFloor = Geo3DUtil.decodeValueFloor(encoded, planetModel);
    assertEquals(minValue, decodedFloor, 0d);
    assertEquals(encoded, planetModel.encodeValue(decodedFloor));
    // Encoding ceiling
    double decodedCeiling = Geo3DUtil.decodeValueCeil(encoded, planetModel);
    assertTrue(decodedCeiling > minValue);
    assertEquals(encoded, planetModel.encodeValue(decodedCeiling));
  }

  public void testMaxValueQuantization() {
    PlanetModel planetModel = randomPlanetModel();
    int encoded = planetModel.MAX_ENCODED_VALUE;
    double maxValue = planetModel.getMaximumMagnitude();
    // Normal encoding
    double decoded = planetModel.decodeValue(encoded);
    assertEquals(maxValue, decoded, 0d);
    assertEquals(encoded, planetModel.encodeValue(decoded));
    // Encoding floor
    double decodedFloor = Geo3DUtil.decodeValueFloor(encoded, planetModel);
    assertTrue(decodedFloor < maxValue);
    assertEquals(encoded, planetModel.encodeValue(decodedFloor));
    // Encoding ceiling
    double decodedCeiling = Geo3DUtil.decodeValueCeil(encoded, planetModel);
    assertEquals(maxValue, decodedCeiling, 0d);
    assertEquals(encoded, planetModel.encodeValue(decodedCeiling));
  }

  // poached from TestGeoEncodingUtils.testLatitudeQuantization:

  /**
   * step through some integers, ensuring they decode to their expected double values. double values
   * start at -planetMax and increase by Geo3DUtil.DECODE for each integer. check edge cases within
   * the double range and random doubles within the range too.
   */
  public void testQuantization() throws Exception {
    Random random = nonAssertingRandom(random());
    PlanetModel planetModel = randomPlanetModel();
    for (int i = 0; i < atLeast(5000); i++) {
      int encoded = random.nextInt();
      if (encoded <= planetModel.MIN_ENCODED_VALUE) {
        continue;
      }
      if (encoded >= planetModel.MAX_ENCODED_VALUE) {
        continue;
      }
      double min = encoded * planetModel.DECODE;
      double decoded = Geo3DUtil.decodeValueFloor(encoded, planetModel);
      // should exactly equal expected value
      assertEquals(min, decoded, 0.0D);
      // should round-trip
      assertEquals(encoded, planetModel.encodeValue(decoded));
      // test within the range
      if (encoded != Integer.MAX_VALUE) {
        // this is the next representable value
        // all double values between [min .. max) should encode to the current integer
        double max = min + planetModel.DECODE;
        assertEquals(max, Geo3DUtil.decodeValueFloor(encoded + 1, planetModel), 0.0D);
        assertEquals(encoded + 1, planetModel.encodeValue(max));

        // first and last doubles in range that will be quantized
        double minEdge = Math.nextUp(min);
        double maxEdge = Math.nextDown(max);
        assertEquals(encoded, planetModel.encodeValue(minEdge));
        assertEquals(encoded, planetModel.encodeValue(maxEdge));

        // check random values within the double range
        long minBits = NumericUtils.doubleToSortableLong(minEdge);
        long maxBits = NumericUtils.doubleToSortableLong(maxEdge);
        for (int j = 0; j < 100; j++) {
          double value =
              NumericUtils.sortableLongToDouble(TestUtil.nextLong(random, minBits, maxBits));
          // round down
          assertEquals(encoded, planetModel.encodeValue(value));
        }
      }
    }
  }

  public void testEncodeDecodeIsStable() throws Exception {
    PlanetModel planetModel = randomPlanetModel();
    int iters = atLeast(1000);
    for (int iter = 0; iter < iters; iter++) {
      double lat = GeoTestUtil.nextLatitude();
      double lon = GeoTestUtil.nextLongitude();

      GeoPoint point = new GeoPoint(planetModel, toRadians(lat), toRadians(lon));

      // encode point
      GeoPoint pointEnc =
          new GeoPoint(
              planetModel.decodeValue(planetModel.encodeValue(point.x)),
              planetModel.decodeValue(planetModel.encodeValue(point.y)),
              planetModel.decodeValue(planetModel.encodeValue(point.z)));

      // encode it again (double encode)
      GeoPoint pointEnc2 =
          new GeoPoint(
              planetModel.decodeValue(planetModel.encodeValue(pointEnc.x)),
              planetModel.decodeValue(planetModel.encodeValue(pointEnc.y)),
              planetModel.decodeValue(planetModel.encodeValue(pointEnc.z)));
      // System.out.println("TEST " + iter + ":\n  point    =" + point + "\n  pointEnc =" + pointEnc
      // + "\n  pointEnc2=" + pointEnc2);

      assertEquals(pointEnc.x, pointEnc2.x, 0.0);
      assertEquals(pointEnc.y, pointEnc2.y, 0.0);
      assertEquals(pointEnc.z, pointEnc2.z, 0.0);
    }
  }

  public void testDistanceQuery() throws Exception {
    Directory dir = newDirectory();
    RandomIndexWriter w = new RandomIndexWriter(random(), dir);

    PlanetModel planetModel = randomPlanetModel();

    // index two points:
    Document doc = new Document();
    doc.add(new Geo3DPoint("field", planetModel, 0.1, 0.1));
    w.addDocument(doc);
    doc = new Document();
    doc.add(new Geo3DPoint("field", planetModel, 90, 1));
    w.addDocument(doc);

    ///// search //////
    IndexReader reader = w.getReader();
    w.close();
    IndexSearcher searcher = newSearcher(reader);

    // simple query, in any planet model one point should be inside and the other outside
    Query q =
        Geo3DPoint.newDistanceQuery("field", planetModel, 0, 0, planetModel.getMeanRadius() * 0.01);
    assertEquals(1, searcher.count(q));

    IOUtils.close(w, reader, dir);
  }

  /**
   * Clips the incoming value to the allowed min/max range before encoding, instead of throwing an
   * exception.
   */
  private static int encodeValueLenient(double x, PlanetModel planetModel) {
    double planetMax = planetModel.getMaximumMagnitude();
    if (x > planetMax) {
      x = planetMax;
    } else if (x < -planetMax) {
      x = -planetMax;
    }
    return planetModel.encodeValue(x);
  }

  private static class ExplainingVisitor implements IntersectVisitor {

    final GeoShape shape;
    final GeoPoint targetDocPoint;
    final GeoPoint scaledDocPoint;
    final IntersectVisitor in;
    final List<Cell> stack = new ArrayList<>();
    private List<Cell> stackToTargetDoc;
    final int targetDocID;
    final int numDims;
    final int bytesPerDim;
    private int targetStackUpto;
    final StringBuilder b;

    // In the first phase, we always return CROSSES to do a full scan of the BKD tree to see which
    // leaf block the document lives in
    boolean firstPhase = true;

    public ExplainingVisitor(
        GeoShape shape,
        GeoPoint targetDocPoint,
        GeoPoint scaledDocPoint,
        IntersectVisitor in,
        int targetDocID,
        int numDims,
        int bytesPerDim,
        StringBuilder b) {
      this.shape = shape;
      this.targetDocPoint = targetDocPoint;
      this.scaledDocPoint = scaledDocPoint;
      this.in = in;
      this.targetDocID = targetDocID;
      this.numDims = numDims;
      this.bytesPerDim = bytesPerDim;
      this.b = b;
    }

    public void startSecondPhase() {
      if (firstPhase == false) {
        throw new IllegalStateException("already started second phase");
      }
      if (stackToTargetDoc == null) {
        b.append("target docID=" + targetDocID + " was never seen in points!\n");
      }
      firstPhase = false;
      stack.clear();
    }

    @Override
    public void visit(int docID) throws IOException {
      assert firstPhase == false;
      if (docID == targetDocID) {
        b.append("leaf visit docID=" + docID + "\n");
      }
    }

    @Override
    public void visit(int docID, byte[] packedValue) throws IOException {
      if (firstPhase) {
        if (docID == targetDocID) {
          assert stackToTargetDoc == null;
          stackToTargetDoc = new ArrayList<>(stack);
          b.append("  full BKD path to target doc:\n");
          for (Cell cell : stack) {
            b.append("    " + cell + "\n");
          }
        }
      } else {
        if (docID == targetDocID) {
          double x = Geo3DPoint.decodeDimension(packedValue, 0, shape.getPlanetModel());
          double y = Geo3DPoint.decodeDimension(packedValue, Integer.BYTES, shape.getPlanetModel());
          double z =
              Geo3DPoint.decodeDimension(packedValue, 2 * Integer.BYTES, shape.getPlanetModel());
          b.append("leaf visit docID=" + docID + " x=" + x + " y=" + y + " z=" + z + "\n");
          in.visit(docID, packedValue);
        }
      }
    }

    @Override
    public Relation compare(byte[] minPackedValue, byte[] maxPackedValue) {
      Cell cell = new Cell(minPackedValue, maxPackedValue);
      // System.out.println("compare: " + cell);

      // TODO: this is a bit hacky, having to reverse-engineer where we are in the BKD tree's
      // recursion ... but it's the lesser evil vs e.g. polluting this visitor API, or
      // implementing this "under the hood" in BKDReader instead?
      if (firstPhase) {

        // Pop stack:
        while (stack.size() > 0 && stack.get(stack.size() - 1).contains(cell) == false) {
          stack.remove(stack.size() - 1);
          // System.out.println("  pop");
        }

        // Push stack:
        stack.add(cell);
        // System.out.println("  push");

        return Relation.CELL_CROSSES_QUERY;
      } else {
        Relation result = in.compare(minPackedValue, maxPackedValue);
        if (targetStackUpto < stackToTargetDoc.size()
            && cell.equals(stackToTargetDoc.get(targetStackUpto))) {
          b.append(
              "  on cell "
                  + stackToTargetDoc.get(targetStackUpto)
                  + ", wrapped visitor returned "
                  + result
                  + "\n");
          targetStackUpto++;
        }
        return result;
      }
    }

    private class Cell {
      private final byte[] minPackedValue;
      private final byte[] maxPackedValue;

      public Cell(byte[] minPackedValue, byte[] maxPackedValue) {
        this.minPackedValue = minPackedValue.clone();
        this.maxPackedValue = maxPackedValue.clone();
      }

      /** Returns true if this cell fully contains the other one */
      public boolean contains(Cell other) {
        for (int dim = 0; dim < numDims; dim++) {
          int offset = bytesPerDim * dim;
          // other.min < this.min?
          if (Arrays.compareUnsigned(
                  other.minPackedValue,
                  offset,
                  offset + bytesPerDim,
                  minPackedValue,
                  offset,
                  offset + bytesPerDim)
              < 0) {
            return false;
          }
          // other.max < this.max?
          if (Arrays.compareUnsigned(
                  other.maxPackedValue,
                  offset,
                  offset + bytesPerDim,
                  maxPackedValue,
                  offset,
                  offset + bytesPerDim)
              > 0) {
            return false;
          }
        }

        return true;
      }

      @Override
      public String toString() {
        double xMin =
            Geo3DUtil.decodeValueFloor(
                NumericUtils.sortableBytesToInt(minPackedValue, 0), shape.getPlanetModel());
        double xMax =
            Geo3DUtil.decodeValueCeil(
                NumericUtils.sortableBytesToInt(maxPackedValue, 0), shape.getPlanetModel());
        double yMin =
            Geo3DUtil.decodeValueFloor(
                NumericUtils.sortableBytesToInt(minPackedValue, 1 * Integer.BYTES),
                shape.getPlanetModel());
        double yMax =
            Geo3DUtil.decodeValueCeil(
                NumericUtils.sortableBytesToInt(maxPackedValue, 1 * Integer.BYTES),
                shape.getPlanetModel());
        double zMin =
            Geo3DUtil.decodeValueFloor(
                NumericUtils.sortableBytesToInt(minPackedValue, 2 * Integer.BYTES),
                shape.getPlanetModel());
        double zMax =
            Geo3DUtil.decodeValueCeil(
                NumericUtils.sortableBytesToInt(maxPackedValue, 2 * Integer.BYTES),
                shape.getPlanetModel());
        final XYZSolid xyzSolid =
            XYZSolidFactory.makeXYZSolid(
                shape.getPlanetModel(), xMin, xMax, yMin, yMax, zMin, zMax);
        final int relationship = xyzSolid.getRelationship(shape);
        final boolean pointWithinCell = xyzSolid.isWithin(targetDocPoint);
        final boolean scaledWithinCell = xyzSolid.isWithin(scaledDocPoint);

        final String relationshipString;
        switch (relationship) {
          case GeoArea.CONTAINS:
            relationshipString = "CONTAINS";
            break;
          case GeoArea.WITHIN:
            relationshipString = "WITHIN";
            break;
          case GeoArea.OVERLAPS:
            relationshipString = "OVERLAPS";
            break;
          case GeoArea.DISJOINT:
            relationshipString = "DISJOINT";
            break;
          default:
            relationshipString = "UNKNOWN";
            break;
        }
        return "Cell(x="
            + xMin
            + " TO "
            + xMax
            + " y="
            + yMin
            + " TO "
            + yMax
            + " z="
            + zMin
            + " TO "
            + zMax
            + "); Shape relationship = "
            + relationshipString
            + "; Quantized point within cell = "
            + pointWithinCell
            + "; Unquantized point within cell = "
            + scaledWithinCell;
      }

      @Override
      public boolean equals(Object other) {
        if (other instanceof Cell == false) {
          return false;
        }

        Cell otherCell = (Cell) other;
        return Arrays.equals(minPackedValue, otherCell.minPackedValue)
            && Arrays.equals(maxPackedValue, otherCell.maxPackedValue);
      }

      @Override
      public int hashCode() {
        return Arrays.hashCode(minPackedValue) + Arrays.hashCode(maxPackedValue);
      }
    }
  }

  public static String explain(
      String fieldName,
      GeoShape shape,
      GeoPoint targetDocPoint,
      GeoPoint scaledDocPoint,
      IndexReader reader,
      int docID)
      throws Exception {

    final XYZBounds bounds = new XYZBounds();
    shape.getBounds(bounds);

    // First find the leaf reader that owns this doc:
    int subIndex = ReaderUtil.subIndex(docID, reader.leaves());
    LeafReader leafReader = reader.leaves().get(subIndex).reader();

    StringBuilder b = new StringBuilder();
    b.append("target is in leaf " + leafReader + " of full reader " + reader + "\n");

    DocIdSetBuilder hits = new DocIdSetBuilder(leafReader.maxDoc());
    ExplainingVisitor visitor =
        new ExplainingVisitor(
            shape,
            targetDocPoint,
            scaledDocPoint,
            new PointInShapeIntersectVisitor(hits, shape, bounds),
            docID - reader.leaves().get(subIndex).docBase,
            3,
            Integer.BYTES,
            b);

    // Do first phase, where we just figure out the "path" that leads to the target docID:
    leafReader.getPointValues(fieldName).intersect(visitor);

    // Do second phase, where we we see how the wrapped visitor responded along that path:
    visitor.startSecondPhase();
    leafReader.getPointValues(fieldName).intersect(visitor);

    return b.toString();
  }
}
